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Abstract 


An  analysis  of  the  three-dimensional  (3-D)  finite  element  formulation  of  Maxwell’s  equations 
governing  classical  electromagnetic  propagation  in  dielectrics  is  given  including  its  analogy  to 
Navier’s  equation.  The  weak  form  of  the  electric  field  equation  is  reviewed  along  with  dispersion 
analysis  and  approximation  equations.  Radiation  boundary  conditions  are  also  explored  to  include 
paraxial  absorber,  Sandler  absorber,  and  other  absorber  comparisons.  In  addition,  time  domains 
vs.  frequency  domains  are  investigated  with  a  listing  of  possible  advantages  and  disadvantages. 
It  was  concluded  that  if  large-scale  calculations  need  to  be  done  today,  time-domain  techniques 
provide  the  most  practicable  means;  however,  it  is  still  premature  to  promote  such  solvers  as 
production  level  tools  for  engineers. 
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1.  Introduction 


The  purpose  of  this  work  is  to  understand  how  the  finite  element  method  or  any  other  numerical 
approximation  may  be  effectively  utilized  in  describing  the  starter  switch  wave  propagation 
phenomena  inherent  in  typical  gas  turbines,  which  are  responsible  for  engine  startup  failures  and 
unnecessary  maintenance  downtimes.  We  will  start  out  by  deriving  and  analyzing  three-dimensional 
(3-D)  transient  finite  element  formulations  of  Maxwell’s  equations  governing  classical 
electromagnetic  propagation  in  dielectrics.  The  derivation  will  provide  a  comprehensive  analytical 
description  of  the  finite  element  equations,  while  the  analysis  will  quantify  accuracy  of  time-harmonic 
plane  wave  propagation  as  a  function  of  wave  direction  and  discretization.  The  derivation  will  be 
limited  to  the  so-called  Cartesian  elements  (hexahedrons)  for  reasons  of  modeling  simplicity  and 
computational  efficiency.  Skewed  elements,  although  an  important  attribute  of  the  finite  element 
method,  will  not  be  discussed  for  the  sake  of  brevity. 

A  tacit  assumption  in  this  formulation  is  that  electromagnetic  waves  are  supported  throughout 
an  infinite  3-D  space.  In  other  words,  there  is  no  outer  boundary  limiting  the  field.  Clearly,  however, 
solving  infinite-domain  field  problems  is  not  practical  Therefore,  3-D  space  must  be  divided  into  a 
finite  interior  or  computational  domain  and  an  infinite  exterior  domain.  By  hypothesis,  the  interior 
includes  all  physical  features  of  interest  while  the  exterior  is  effectively  “featureless”  so  that  it 
propagates  energy  outward  with  negligible  backscatter.  This  artifice  allows  replacement  of  the 
exterior  domain  with  a  so-called  radiation  boundary  condition  on  the  outer  surface  of  the  interior 
domain  and  defines  the  boundary  value  part  of  the  problem. 

The  finite  element  formulation  of  initial-boundary  value  problems  governed  by  time-dependent 
partial  differential  equations  (PDEs) — Maxwell’s  equations  in  particular— consists  of  the  following 
formal  steps: 

1.  Partition  the  problem’s  interior  domain  into  a  number  of  logically  regular,  contiguous 
subdomains. 
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2.  Represent  the  field  over  each  subdomain  by  a  simplified  basis  function  that  interpolates 
between  discrete  field  points  or  nodes. 

3 .  Convert  the  pointwise  partial  differential  operator  to  an  equivalent  but  “weaker”  scalar  integral 
operator  admitting  lower  order  derivatives. 

4.  Evaluate  the  integral  operator  for  the  simplified  field  basis,  giving  an  algebraic  system  of 
equations  on  the  nodal  field  vector  and  its  time  derivatives. 

5.  Apply  a  radiation  condition  on  the  interior  domain’s  boundary  in  order  to  simulate  scattering 
into  the  infinite  exterior  domain. 

6.  Solve  the  system  of  ordinary  differential  equations  (ODEs)  in  time  using  finite  differences, 
modal  analysis,  etc. 

The  finite  element  part,  steps  1-4,  yields  an  approximate  integration  of  the  PDE’s  spatial 
differential  operator.  In  other  words,  the  so-called  finite  element  discretization  is  nothing  more  than 
a  quadrature  formula.  The  remaining  one-dimensional  (1-D)  temporal  problem,  step  6,  is  typically 
integrated  in  a  more  conventional  fashion. 

Formal  reduction  of  the  pointwise  partial  differential  equation  to  a  finite  element  form  may  be 
accomplished  in  at  least  two  ways — the  method  of  weighted  residuals  (Galerkin’s  method)  and  a 
variational  principle.  They  are  fundamentally  equivalent  for  self-adjoint  (symmetric)  differential 
operators,  provided  that  consistent  assumptions  are  made,  although  the  method  of  weighted  residuals 
is  more  general  The  approach  chosen  usually  depends  on  the  analyst’s  perspective  or  background, 
e.g.,  mathematical,  engineering,  etc.,  rather  than  any  compelling  analytical  reason.  Given  here  is  a 
somewhat  simplified  derivation  that  is  comprehensive  yet  minimizes  nomenclature  and  historical 
biases.  For  further  reading,  a  succinct  description  of  the  formalism  in  the  context  of  linear  differential 
operators  may  be  found  in  Zienkiewicz  (1977). 
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The  most  difficult  part  of  the  finite  element  formulation  of  propagation-type  problems  is  deriving 
an  effective  radiation  or  absorbing  boundary  condition.  This  is  an  approximate  condition  on  the 
exterior  boundary  of  the  finite  element  model  that  discriminates  between  incident  (illumination)  and 
scattered  radiation  and  selectively  absorbs  the  scattered  part,  mimicking  radiation  into  an  infinite, 
nonreflecting  exterior  domain.  An  effective  condition  that  is  sufficient  for  simultaneous  plane  wave 
illumination  and  scattering  will  be  attempted  in  the  succeeding  parts  of  this  treatise. 

2.  Maxwell’s  Equations 

Maxwell’s  equations  provide  the  mathematical  basis  for  rigorous  analysis  of  classical 
electromagnetic  wave  propagation.  In  particular,  they  provide  a  complete  description  of  macroscopic 
optical  phenomena  in  dielectric  media.  There  are  lower  limits  on  size  and  intensity  where  quantum 
behavior  becomes  significant,  but  for  nearly  all  scales  of  practical  interest,  Maxwell’s  equations  are 
both  necessary  and  sufficient.  An  attempt  at  the  various  forms  of  the  equations  along  with  ancillary 
relations  that  are  useful  for  numerical  algorithm  development  will  be  made  in  the  most  terse  and 
sufficient  manner. 

Maxwell  originally  proposed  an  arcane  system  of  20  equations  in  20  unknowns.  The  system 
was  subsequently  simplified  by  Heaviside  and  Hertz  to  its  modem  form,  namely 

-B=VxE,  D  +/  =  Vxff,  (1) 

where  B  is  magnetic  induction,  E  is  electric  field  intensity,  D  is  electric  displacement,  J  is  current 
density,  and  H  is  magnetic  field  intensity.  Continuity  of  E  and  H  is  required  in  order  to  define  the 
spatial  derivatives.  Note  that  in  (1)  bold  letters  represent  vectors,  Vx  is  the  curl  operator,  and  time 
derivatives  are  denoted  by  a  dot  above  the  variable. 
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A  useful  alternative  to  Maxwell’s  pointwise  PDEs  are  volumetric  forms  obtained  by  integrating 
(1)  over  space.  Applying  the  so-called  curl  theorem — analogous  to  the  divergence  theorem  of 
Gauss — to  the  resulting  integrals  of  V  x  E  and  Vxff  yields  the  vector  integral  equations 

-fBdQ  =  f  n  x  EdL  ,  J(D  +  J)dQ  =  j  n  x  HdZ  ,  (2) 

Q  S  0  S 

where  Q  is  the  domain  of  integration,  S  is  its  boundary,  and  n  is  the  outward  unit  normal  vector  to 
2.  This  form  removes  the  restriction  on  field  continuity  since  spatial  derivatives  no  longer  appear. 
More  conventional  scalar  integral  equations  may  be  written  by  applying  Stokes’  theorem  to 
Maxwell’s  equations,  yielding  the  famous  laws  of  Ampere  and  Faraday. 

To  make  Maxwell’s  equations  determinate  for  B,  E,  D,  J,  and  H,  so-called  constitutive  relations 
must  be  defined.  In  nearly  all  cases,  the  linear  relations 

B  =  ptf,  D  =  zE,  J  =  oE  (3) 

suffice.  Proportionality  factors,  p,  e,  and  o,  are  magnetic  permeability,  dielectric  permittivity,  and 
conductivity,  respectively.  Provided  the  medium  is  isotropic,  these  factors  are  scalars,  otherwise  they 
are  tensors.  Substituting  (3)  into  (1)  gives  the  determinate  form  of  Maxwell’s  PDEs, 

-ptf  =  V  x  E,  eE  +  oE  =  V  x  H,  (4) 

relating  the  time  rate  of  change  of  the  magnetic  field  to  the  curl  (vorticity  or  “swirl”)  of  the  electric 
field  and  vice  versa.  These  assume  that  e  and  p  vary  with  time  slowly,  if  at  all,  in  comparison  to  the 
fields  themselves.  For  nonmagnetic  materials,  permeability  p  is  essentially  equal  to  its  vacuum 
value,  p0,  everywhere. 

Since  E  is  the  primary  field  unknown  for  analyses  in  dielectric  media,  H  may  be  eliminated 
between  the  two  curl  equations  in  (4)  and  treated  as  a  secondary  or  derived  quantity.  Taking  the  time 
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derivative  of  the  second  equation  in  (4),  the  curl  of  the  first,  and  eliminating  the  term  with  H  gives 
the  second  order  partial  deferential  equation, 


eE  +  —  V  xV  x  E,  (5) 

P 

where  the  constant  magnetic  permeability  is  brought  outside  the  curl  operator.  Such  a  simple 
equation  cannot  be  written  for  H  because  e,  unlike  p,  is  a  function  of  position  and  its  gradient  must 
be  included.  Using  the  vector  identity,  VxVx£  =  V(V>£)-V2£,  equation  (5)  may  also  be 
written 


eE  +  o£=-^E-  -  V(V-£). 
P  P 


Note  the  similarity  between  this  equation  and  Navier’s  equation 


(6) 


pi)  =  C&U  -  (A  +  G)  V  (V  •  17) 


(7) 


describing  displacement  U  in  a  linear  elastic  medium. 

It  should  be  noted  for  completeness  that  the  vector  fields  in  (1)  are  ultimately  caused  by  some 
distribution  of  electric  charge,  p,  and  current,  J,  generally  related  by  the  continuity  equation, 

V'7  +  p=  0,  (8) 

expressing  pointwise  conservation  of  charge.  Taking  the  divergence  of  (1),  substituting  (8),  and 
integrating  over  time,  assuming  a  quiescent  initial  or  final  state,  gives  the  divergence  conditions, 

V  B  =  0,  V-Z)  =  p.  (9) 

These  are  often  appended  to  Maxwell’s  equations  but  are  dependent  conditions,  either  in  part 
(V  •  B  =  0),  or  wholly  if  charge  is  conserved  (V  •  D  =  p).  Observe  that  by  solving  V  •  D  =  V  •  zE  =  p 
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for  V  •  E  =  p/e  -  V  e/e  •  E  and  replacing  V  •  E  in  (6),  a  second  order  equation  results  that  implicitly 
incorporates  charge  conservation. 


3.  A  Weak  Form  of  the  Electric  Field  Equation 


To  apply  the  conventional  finite  element  formalism  to  Maxwell’s  equations,  it  is  convenient  to 
start  with  the  second  order  PDE  on  the  electric  field,  (5),  rather  than  the  original  system  of  first  order 
equations,  (4).  Strict  solutions  of  this  equation  must  possess  at  least  second  derivatives;  however, 
it  is  impractical  to  require  such  continuity  from  numerical  approximations.  A  better  approach  is  to 
rewrite  the  equation  in  an  integral  form  admitting  lower  order  derivatives.  This  is  the  so-called  weak 
formulation. 

To  derive  the  weak  form  of  (5)  it  is  necessary  to  define  another  field  over  the  wave  domain,  the 
so-called  test  function,  G(x,t).  This  is  a  completely  arbitrary  function  within  wave  domain  Q.  Taking 
the  inner  (dot)  product  of  (5)  with  G  and  integrating  over  Q  gives 

fG-(eE  +  oE)dQ  =  -Jg- —VxVxEdQ.  (10) 

Multiplication  by  a  test  function  and  integration  reduces  the  pointwise  vector  equation  to  a  volumetric 
scalar  equation — the  weak  form.  It  is  easy  to  prove  the  assertion  that  if  this  integral  equation  is 
satisfied  for  any  G,  then  the  PDE  is  necessarily  satisfied  at  all  points  in  the  domain.  The  converse  is 
certainly  true;  but  if  the  PDE  is  not  satisfied  in  some  subdomain,  then  a  test  function  can  be  chosen 
that  makes  the  integral  nonzero,  hence  the  assertion  is  true. 

Consider  the  right-hand  integrand  in  (10),  i.e.,  G  •  V  x  (V  x  E)  after  factoring.  From  the  vector 
identity,  V  •  (A  x  B)  =  B  *  (V  x  A)  -  A  •  (V  x  B),  this  integrand  can  be  written 

G*V  x  V  x  E  =  V  x  E-7  x  G  -  V*G  xVx£.  (11) 
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Integrating  and  applying  the  divergence  theorem  to  the  second  term  gives 


jG-VxVxEdQ  =  J  V  x  2?  •  V  x  GdQ  +  f  G-  nxVx  Edh .  (12) 

o  os 

In  the  surface  integral,  n  is  the  outward  unit  normal  to  2  and  the  integrand  has  been  rearranged 
according  to  the  rule  for  scalar  triple  products.  This  identity  is  the  vector  analog  of  Green’s  identity, 
e.g.,  see  Stratton  (1941),  and  is  simply  the  result  of  multidimensional  integration  by  parts. 
Substituting  (3)  into  (3.1),  the  volume-averaged  scalar  equation  becomes 

jG-(eE  +  a£)dQ  =  - -J"VxE’VxGdQ- - jG-nxVxEdS .  (13) 

o  ^  q  ^  s 

The  critical  result  expressed  in  (13)  is  that  the  volume  integral  of  the  second  order  spatial  operator 
has  been  replaced  by  “weaker”  volume  and  surface  integrals  of  first  order  operators. 


4.  Reduction  to  an  Ordinary  Differential  Equation 


The  basis  for  transforming  the  volumetric  partial  differential  equation  to  an  ordinary  differential 
equation  is  an  assumption  on  the  mathematical  form  of  wave  fields  in  domain  Q.  In  particular, 
fields  are  assumed  to  be  separable  in  space  and  time,  namely, 

E(x,t)  =S(x)f(t),  G(x,t)  =S(x)g(t),  (14) 

where  x  is  the  position  vector,  matrix  S(jc)  represents  the  field’s  spatial  variation,  and  vector /(t)  or 
g( t)  represents  the  time  variation.  Note  that  the  same  spatial  variation  is  assumed  for  E  and  G  in  (14). 
This  assumption,  associated  with  the  name  of  Galerkin  in  the  finite  element  literature,  is  particularly 
convenient  because  it  yields  a  symmetric  system  of  equations. 

Separable  representation  (14)  may  be  interpreted  in  a  number  of  ways.  For  example,  S(x)  can  be 
eigenvectors  (mode  shapes),  whence /(f)  are  the  corresponding  eigenvalues  (frequencies)  for  the 
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given  domain;  or  S(x)  can  be  a  multidimensional  Fourier  series,  for  which/(r)  are  the  time-dependent 
Fourier  coefficients.  Alternatively,  it  is  easy  to  show  that  (14)  is  the  functional  form  of  a 
multidimensional  Taylor  series.  In  that  case,  fit)  is  an  infinite  vector  of  all  derivatives  at  point  x  and 
S(x)  is  a  corresponding  infinite  matrix  consisting  of  the  coefficients  of  these  derivatives  in  the  series, 
i.e.,  powers  of  the  local  space  coordinates.  Note  that  each  of  these  interpretations  can  provide  a  local 
or  global  basis  for  analysis. 

Substituting  separable  solutions  (14)  into  the  integrands  in  (13)  and  applying  the  vector 
differential  operator  to  those  on  the  right  side  gives 

G‘tE  =  gTSTeSf  ,  G  o£=grSToSf 
Vx£ *VxG  =  £r(Vx  S)r(VxS)/ 

G  nx  (yx.E)=gTSrn  x (Vx 5)/.  (15) 

Therefore,  substituting  into  (13),  moving  the  vector  functions  of  time  outside  the  integrals,  and 
rearranging  yields 


gT\Mf+Cf-Kf-Bj)  =  0, 


(16) 


where 


M 


=  fSTeSdQ,  C-  JS  ToSdQ,  K  =  ~j  (VxS)r(VxS)dQ 


(17) 


are  symmetric  coefficient  matrices  defined  by  the  volume  integrals  and 


B  = 


-1  jSTnxVxSd2 


(18) 


is  the  matrix  defined  by  the  surface  integral.  Since  g(t)  is  arbitrary,  (16)  is  satisfied  when 


S 


Mf  +  Cf  ={K  +  B)f. 


(19) 


This  is  the  global  ordinary  differential  equation  equivalent  to  Maxwell’s  PDEs  in  Q.  Of  course,  the 
utility  of  (19)  depends  on  the  choice  of  separable  field  representation,  i.e.,  S(x)  and  fit). 

5.  The  Finite  Element  Equations 


Given  the  previous  mathematical  preamble,  the  finite  element  procedure  consists  of  partitioning 
or  discretizing  interior  domain  into  a  number  of  subdomains  or  finite  elements.  The  field  is 
approximated  over  each  element  by  an  interpolating  or  shape  function  depending  on  values  at  discrete 
nodes  on  or  in  the  element.  This  yields  a  convenient  local  basis  (in  contrast  to  a  global  basis)  for 
evaluating  the  model’s  matrix  coefficients  in  (17)  using  an  element-by-element  summation. 

To  provide  some  degree  of  field  continuity  across  element  boundaries,  most  of  the  discrete  nodes 
are  defined  on  the  element  surfaces  and  shared  by  adjacent  elements.  Provided  that  they  cover  the 
domain,  the  elements  and  shape  functions  may  be  completely  arbitrary.  However,  it  is  advantageous 
in  terms  of  modeling  and  computation  to  make  elements  as  simple  as  possible.  Thus,  in  three 
dimensions,  simple  hexahedron  or  brick  shapes,  i.e.,  so-called  Cartesian  elements,  are  favored  with 
low-degree  shape  functions  based  on  nodes  at  the  comers,  and  on  the  faces  and  edges  as  the  degree 
of  interpolant  is  increased.  The  tri-linear  shape  function  (linear  in  each  direction)  with  comer  nodes 
is  the  lowest  degree  that  provides  field  continuity  in  3-D;  hence,  this  is  the  most  “elemental” 
interpolant.  A  simple  analog  to  this  shape  function  is  the  common  trapezoidal  rule  used  in  numerical 
integration  of  functions  in  one  or  more  dimensions. 

With  the  domain  covered  by  an  assemblage  of  elements  and  nodes,  individual  elements,  m  =  1, 
M,  and  nodes,  n  =  1,  N,  are  consecutively  numbered  in  some  convenient  fashion.  The  global 
unknown  vector,  f(t)  in  (14),  is  written  as  f(t)  =  [/i,/2,/3]T,  where  vectors/^  k=  1,2, 3  are  the  three 
components  of  the  electric  field  at  the  N  ordered  nodes  of  the  assemblage. 
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We  consider  a  single,  eight-node,  Cartesian  finite  element  or  brick  with  element  numbering  and 
coordinates  system.  This  is  the  generic  element  composing  any  assemblage,  with  element  node 
numbers  related  to  global  node  numbers  by  a  simple  map.  The  shape  function  matrix  and  node  vector 
for  element  m  are  written  as 


Sm(x) 


sOT(x)  0  0 

I 

0  sm(x )  0 

0  0  sm(x) 


fm(t) 


(20) 


where  sm(x)  and  0  are  row  8-vectors  and  fkm(t)  are  column  8-vectors  for  the  three  field  components. 
The  same  spatial  basis,  Le.,  sm  (x),  is  assumed  for  each  field  component.  Note  that  the  curl  term  in 
(13)  becomes 


VxSm(x) 


0 

dsm 

dz 


dsm  dsm 
dz  dy 


0 


dsm 

dx 


dsm  dsm 
dy  dx 


The  canonical  tri-linear  shape  function  assumed  here  is  the  row  vector 


s 


m  _ 


i 

8 


’  ( 1  -  2x  /  Ax)  ( 1  -  2y/ Ay )  ( 1  -  2z  /  Az) ' 
(l+2x/Ajt)(l-2y/Ay)(l-2z/Az) 
(1  +2x  /  Ax)(l  +2y  /  Ay)(l-2z  /  A  z) 
(l-2x/  Ax)(l  +2y  /  Ay)(l-2z  /  A  z) 
(1  -  2x  /  Ax)(  1  -  2y  /  Ay)(  1 +2z  /  Az) 
(1  +2x  /  Ax)(l-2y  /  Ay)(l  +2z  /  A  z) 
(1  +2x  /  Ax)(l  +2y  /  Ay)(l  +2z  /  A  z) 
(l-2x  /  Ax)(l  +2y  /  Ay)(l+2z  /  A  z)_ 


r 


(21) 


(22) 
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Each  component  (1  to  8)  of  this  vector  is  unity  at  the  corresponding  node  (1  to  8)  and  decreases 
linearly  to  zero  at  all  other  nodes. 

Given  definite  forms  of  the  element  shape  function  and  nodal  vector  /,  the  elemental  matrix 
coefficients  can  be  computed.  These  are  found  by  substituting  the  above  into  integral  definitions  (17) 
and  evaluating.  For  most  modeling  purposes,  e  and  o  may  be  assumed  constant  over  an  element; 
hence,  they  are  typically  factored  out  of  the  integration.  There  is  no  need  to  evaluate  surface  integral 
term  Bm  since  it  is  irrelevant  at  the  element  level,  although,  as  a  point  of  interest,  it  is  numerically 
equal  to  -K"1  by  virtue  of  (1 1)  and  (13). 

The  algebra  and  integration  of  the  24  x  24  matrices  are  somewhat  tedious  and  best  done 
symbolically  using  a  program  like  Macsyma  or  Mathematica.  All  of  the  evaluations  described  in 
the  following  sections  were  done  in  Mathematica  (Wolfram  1988)  because  of  its  general  utility  and 
availability  on  personal  computers. 

Rather  than  listing  the  24  x  24  coefficient  matrices,  suffice  it  to  say  that  they  are  fully  populated, 
symmetric,  and  real,  with  Mm/e  =  Cmo  when  e  and  o  are  constant.  With  the  element  coefficient 
matrices  thus  evaluated,  the  global  system  of  finite  element  equations  is  assembled  by  inflating  the 
element  matrices,  i.e.,  mapping  the  element  row  and  column  positions  to  the  global  row  and  column 
positions  and  summing  the  contribution  from  each  element  in  the  assemblage. 

6.  Dispersion  Analysis  of  the  Finite  Element  Equations 

The  numerical  solution  of  PDEs  is  not  as  straightforward  as  the  previous  derivations  would 
suggest.  There  are  many  variations  on  the  basic  numerical  theme  that  need  to  be  considered  in  the 
context  of  the  physical  application.  An  effective  means  of  assessing  pros  and  cons  is  dispersion 
analysis.  It  is  particularly  useful  for  quantifying  wave-analytic  properties  of  the  numerical 
solutions. 
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Since  any  finite  element  discretization  introduces  an  artificial  length  scale  (element  size),  wave 
propagation  through  the  model  is  necessarily  dispersive,  i.e.,  phase  velocity  depends  on  frequency. 
Also,  properties  exhibit  directional  dependence  by  virtue  of  the  element’s  shape,  so  phase  velocity 
is  also  anisotropic.  This  dispersive  and  anisotropic  behavior  is  unavoidable  in  any  discrete  numerical 
solution.  Fortunately,  errors  can  be  made  negligible  for  practical  purposes  by  ensuring  that  element 
size  is  small  compared  to  the  minimum  wavelength  in  the  propagating  signal. 

A  simplified  model  to  quantify  these  errors  consists  of  the  2  x  2  x  2  assemblage  or  molecule  of 
Cartesian  elements  for  the  case  of  square  elements.  Since  nodes  are  only  coupled  to  their  nearest 
neighbor,  this  simple  molecule  is  sufficient  to  write  the  complete  set  of  equations  governing  the 
electric  field  at  the  interior  node.  These  equations  provide  the  dispersion  relations  that  completely 
describe  wave-analytic  properties  of  an  infinite  Cartesian  grid. 

To  assemble  the  molecule,  a  left-to-right,  top-to-bottom,  front-to-back  numbering  convention  is 
used.  The  three  layers  or  sheets  of  nodes  and  the  two  sheets  of  elements  are  so  numbered,  starting 
from  the  upper  left  comer.  The  node  numbering  sequence  for  a  single  element  is  mapped  to  the 
molecule  node  numbering  for  each  element  using  a  simple  lookup  table.  Inflating  the  element 
equations  via  this  mapping  and  assembling  them  by  summation  yields  81  x  81  M  and  K  matrices  in 
the  finite  element  model,  Mf  -  Kf  =  0.  Unknown  vector  f(t)  is  composed  of  27  x-field,  27  y-field, 
and  27  z-field  unknowns  at  the  27  nodes  of  the  molecule.  Equations  for  the  center  node  may  be  given 
by  rows  14, 41,  and  68. 

We  assume  that  a  time-harmonic,  linearly  polarized,  plane  wave  propagates  through  the  model 
in  the  direction  of  unit  wave  vector  k  =  {sin  0  cos  <j>,  sin  0  sin  <f>,  cos  ©}T,  specified  by  spherical 
angles,  0  and  4>.  This  is  given  by 


E  =  Asin{T3(k  •  x/  v- 1)} ,  (23) 

where  A  is  the  unit  polarization  vector,  x  is  the  position  vector,  and  v  is  the  phase  velocity.  Given  this 
prescription  of  the  incident  field  vector, /(f)  is  found  directly  by  evaluating  (23)  at  the  x  coordinate 
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of  each  of  the  27  nodes.  Instead  of  simply  differentiating  (23)  to  obtain  /,  it  is  evaluated  from  the 
central  difference  approximation 

E(i)  *  (E(t  +  At)  -  2E(t)  +  E(t  -  At))  /  At2 ,  (24) 

corresponding  to  the  typical  stepwise  forward  integration  scheme  with  timestep  At  and  error  on  the 
order  of  At4. 

The  problem  is  to  determine  phase  velocity  v  compatible  with  the  finite  element  equations  and 
assumed  form  of  the  propagating  field.  This  is  accomplished  by  evaluating  solution  vectors/and/ 
at  the  nodes  from  (23)  and  (24),  isolating  the  three  field  equations  for  the  interior  node  in  M/\  Le., 
rows  14,  41,  and  68,  and  solving  for  v  and  A.  It  is  convenient  to  rewrite  the  three  equations  as  VA 
=  0  where  V  is  a  3  x  3  matrix  coefficient.  This  system  of  equations  has  a  solution  if  and  only  if  the 
determinant  of  V  vanishes,  i.e.,  W I  =  0,  yielding  a  nonlinear  scalar  equation  on  phase  velocity  v.  This 
is  the  dispersion  relation  and  its  solution  near  c  =  1  /  /ep  must  be  found  numerically.  Substituting 
v  back  into  V  gives  the  homogeneous  system  of  equations  governing  polarization  vector  A.  In 
particular,  A  is  the  space  mapped  by  linear  transformation  V  into  the  null  vector,  i.e.,  A  is  equal  to 
the  null  space  of  V.  These  solutions  are  described  below. 

Near  v  =  c  the  dispersion  relation  is  found  to  be  second  order  in  general,  i.e.,  it  looks  like  a 
parabola  locally.  Therefore  it  exhibits  two  solutions  near  c.  This  multiplicity  can  be  traced  to  the 
element  shape  function  row  vector,  equation  (22).  Evaluating  the  products  yields  terms  proportional 
to  1,  x,  y,  z,  xy,  xz,  yz,  xyz  for  each  vector  component,  corresponding  to  the  first  eight  terms  in  a 
Taylor  series  (recall,  there  are  eight  nodal  values;  hence,  eight  terms  in  the  series  are  determinate). 
However,  by  a  simple  rotation  of  coordinates  ( X ,  Y,  Z  say),  the  last  four  terms  may  be  converted  to 
product  terms  like  X2,  Y2,  Z2,  and  X3,  Y3,  Z*.  Therefore,  the  shape  function’s  form  is  not  rotationally 
invariant,  unlike  the  original  differential  equation.  This  “rotational  variability”  shows  itself  as  two 
phase  velocities  in  each  direction. 
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Numerical  dispersion  results  for  the  case  of  a  cubic  Cartesian  grid,  showing  v/c  vs.  spherical 
angles  of  incidence  and  maximum  v/c  vs.  discretization.  Discretization  is  measured  by  the  number 
of  elements  supporting  a  wavelength.  Shape  of  the  normalized  v  surface  is  nearly  independent  of 
discretization.  Note  that  maximum  error  occurs  along  the  Cartesian  axes  and  phase  velocity  is 
always  greater  than  c.  See  Sandler  (1991)  for  a  more  complete  discussion  and  explanation  of  the 
aforementioned  results. 

Eigenanalysis  of  V  shows  that,  in  general,  this  3x3  matrix  is  doubly  degenerate,  i.e.,  there  is  one 
nonzero  eigenvalue  and  two  “zero”  eigenvalues  that  are  clustered  very  near  zero.  Eigenvectors  of 
these  “zero”  eigenvalues  span  the  nullspace  of  V — a  plane  in  this  case.  Vector  A  lies  in  this  so-called 
polarization  plane,  which  is  perpendicular  to  the  eigenvector  of  the  nonzero  eigenvalue  of  V, 
denoted  q.  Since  Maxwell’s  equations  represent  transverse  waves,  ideally  the  polarization  plane 
is  perpendicular  to  wave  vector  k,  i.e.,  q  and  k  are  colinear;  however,  grid  anisotropy  produces  a 
“transversality”  error,  which  can  be  measured  by  the  angle  between  q  and  k. 

7.  Approximate  Finite  Element  Equations 

Dispersion  analysis  of  the  exact  finite  element  equations  indicates  four  difficulties.  The  first  is 
philosophical,  namely,  that  the  numerical  phase  velocity  in  a  vacuum  is  greater  than  the  continuum 
speed  of  light — which  physical  objects  can  never  exceed.  It  would  be  preferable  to  have  a 
conservative  solution  where  numerical  velocities  are  always  less  than  continuum  velocities. 

The  second  difficulty  is  that  phase  velocity  errors  are  greatest  in  the  local  coordinate  directions, 
i.e.,  normal  to  the  element  faces.  In  practice,  it  is  natural  to  align  numerical  models  with  these 
coordinates  and  also  to  gather  information  on  waves  traveling  along  or  near  them,  e.g.,  for  imaging. 
Therefore,  it  would  be  better  if  propagation  errors  were  minimized  rather  than  maximized  in  these 
directions. 
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The  third  difficulty  is  existence  of  two  waves  in  any  direction.  This  is  typical  of  anisotropic 
media,  where  the  two  waves  are  denoted  ordinary  and  extraordinary  waves.  Since  the  phase 
velocities  are  close,  the  two  waves  are  numerically  indistinguishable  in  most  cases.  Nonetheless,  from 
an  analytical  viewpoint,  particularly  with  respect  to  boundary  conditions,  this  is  an  unwelcome 
complication. 

The  fourth  difficulty  is  the  number  of  floating  point  operations  necessary  to  evaluate  and  solve 
the  exact  finite  element  equations.  Acording  to  Sandler  (1991),  there  is  at  least  a  factor  of  5  more 
than  necessary  for  a  conventional  finite  difference  approximation  of  Maxwell’s  equations.  This 
profligacy  is  a  deterrent  to  finite  elements  despite  their  marked  advantage  for  modeling  geometrically 
complicated  features. 

These  difficulties  are  known  to  numerical  wave  propagation  analysts,  particularly  in  the  elasticity 
community.  One  approach  that  claims  to  remedy  all  of  them  is  approximation  of  the  coefficient 
matrices  by  parameter  lumping  and  reduced  integration.  Parameter  lumping  diagonalize  the  M  and 
C  matrices  by  placing  the  sum  of  each  row  in  the  diagonal  position  and  zeroing  the  off-diagonal 
terms.  Reduced  integration  applies  to  evaluation  of  the  K  matrix,  using  an  approximate  quadrature 
rule,  e.g.,  the  simple  rectangular  rule,  or  more  generally,  single-point  Gaussian  quadrature. 

Dispersion  analysis  of  the  lumped  parameter,  reduced  integration  finite  element  equations  is  done 
in  the  same  way  as  for  the  exact  equations.  The  dispersion  relation  is  still  locally  parabolic  but  just 
reaches  zero  at  its  maximum  (or  minimum  depending  on  the  sign  chosen)  so  there  is  only  one  phase 
velocity.  Velocity  dependence  on  angle  of  incidence  is  qualitatively  similar,  but  numerical  phase 
velocities  are  always  less  than  exact  phase  velocities.  Maximum  error  occurs  near  the  space  diagonal 
(0  =  45°,  <f>  -  45°),  and  numerous  authors  show  that  absolute  phase  velocity  error  is  only  slightly 
greater  than  that  exhibited  by  the  exact  equations. 

Clearly,  the  approximate  finite  element  equations  solve  the  first  three  problems  mentioned  at  the 
beginning  of  this  section.  A  count  of  arithmetic  operations  shows  that  they  also  require  less  than 
half  the  operations  needed  to  evaluate  the  exact  equations.  The  physical  basis  for  these 
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approximations  can  be  found  in  the  differential  and  integral  forms  of  Maxwell’s  equations,  (2).  In 
particular,  the  lumped  parameter  and  reduced  integration  finite  element  equations  can  be  shown  to 
follow  from  the  first  of  (1)  and  the  second  of  (2)  written  as 


H  =  -  —VxE  ,  I eMl  =  I nxHdZ 


for  the  case  of  zero  conductivity. 


(25) 


8.  Radiation  Boundary  Conditions 


Discrete  numerical  methods  like  finite  elements  or  finite  differences  are  necessarily  formulated  on 
finite  spatial  grids — whether  the  actual  domain  being  modeled  is  finite  or  infinite.  This  domain 
truncation  introduces  artificial  boundaries  that  must  be  treated  with  special  care  in  order  to  minimize 
nonphysical  wave  reflections.  These  trap  energy  that  would  otherwise  be  radiated  and  establish 
undesirable  resonances  within  the  grid.  Note  that  this  is  true  regardless  of  the  solution  scheme 
applied,  in  either  the  time-  or  frequency-domain.  The  solution  to  the  problem  is  to  apply  so-called 
radiation  or  absorbing  boundary  conditions  on  the  model’s  exterior  surfaces. 

In  optics-type  problems,  there  is  an  added  complication  because  some  form  of  illumination  is 
usually  prescribed  over  the  model.  This  is  accommodated  by  an  illumination  boundary  condition, 
designed  to  apply  the  known  incident  electromagnetic  field  on  the  model  boundaries  and  transmit  or 
absorb  the  scattered  field  as  if  the  model  extended  out  to  infinity.  Determination  of  an  incident 
electric  field  consistent  with  the  numerical  propagation  characteristics  of  the  grid  is  a  nontrivial 
calculation,  especially  in  the  presence  of  arbitrary  angles  of  incidence  and  model  topography.  These 
numerical  problems  are  addressed  elsewhere  in  this  report.  Here,  it  is  assumed  that  the  incident  field 
is  known,  and  approximate  radiation  conditions  for  the  scattered  fields  are  described. 

In  formulating  radiation  boundary  conditions,  there  is  a  tradeoff  between  accuracy  and 
complexity.  Generally,  the  more  boundary  nodes  that  are  coupled  by  the  boundary  formulation,  the 
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more  accurate  and  computationally  intensive  the  condition  will  be  for  a  given  model  size. 
Increasingly  accurate  radiation  conditions  can,  in  principle,  allow  smaller  and  smaller  models  around 
the  scattering  features  of  interest.  Development  along  these  lines  may  well  be  warranted,  but  for  this 
report,  attention  is  restricted  to  absorbing  conditions  that  are  local  in  space  and  time  for  compatibility 
with  the  explicit  time  integration  approach.  In  this  way,  no  more  than  a  small  fraction  of  the 
computational  effort  is  expended  on  the  evaluation  of  boundary  conditions. 

Here,  the  basis  for  radiation  conditions  is  the  paraxial  wave  equation,  i.e.,  an  equation  valid  for 
propagation  in  (and  around)  a  selected  direction.  The  prototypical  example  is  the  1  -D  wave  operator, 
which  can  be  factored  as 


a2 

2  a2  ' 

-  rL  = 
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j  +  /*  _ 
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L  -  li 

dx\ 

jar  dx] 
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where  the  factors 

.i.±ci.\E  =  °  (27) 

dr  dx] 

represent  the  right  and  left  traveling  waves.  Thus,  if  the  time  and  space  derivatives  are  related 
according  to  this  operator,  at  the  ends  of  a  1-D  domain,  then  an  “exact”  radiation  boundary  condition 
results;  it  cannot  be  truly  exact  in  practice  because  approximations  are  implicit  in  the  numerical 
implementation.  This  is  the  well-known  normal  incidence  condition,  the  so-called  “black”  boundary. 
It  turns  out  to  be  quite  reasonable  for  many  applications,  particularly  for  boundaries  very  far  from  a 
compact  body  in  a  homogeneous  material,  e.g.,  for  conventional  radar-like  scattering  problems. 

In  order  to  significantly  enhance  boundary  condition  performance,  higher  order  paraxial 
approximations  of  the  multidimensional  Maxwell’s  equations  must  be  derived  and  applied  as 
discussed  below.  Also,  a  condition  proposed  by  Sandler  (1991)  is  described.  The  4th  order  paraxial 
and  Sandler  conditions  are  demonstrated  to  be  roughly  equivalent  in  accuracy,  although  certain 
implementation  issues  appear  to  favor  the  latter. 
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In  their  paper,  “Absorbing  Boundary  Conditions  for  Acoustic  and  Elastic  Wave  Equations,” 
Clayton  and  Engquist  (1977)  developed  an  absorbing  boundary  condition  for  the  two-dimensional 
(2-D)  elastic  wave  equation  using  a  paraxial  approximation  method.  This  approach  has  been 
extended  here  to  Maxwell’s  equations  with  zero  conductivity.  For  the  sake  of  brevity,  only  the 
resulting  paraxial  equations  are  stated,  since  the  derivation  is  not  immediately  relevant.  Engquist 
[1991]  claims  that  extension  to  the  case  of  finite  conductivity  should  also  be  possible. 

The  paraxial  versions  of  Maxwell’s  equations  are  derived  from  the  equivalent  second  order  partial 
differential  equation  on  E,  from  equation  (5).  In  the  following  equations,  it  is  assumed  that  the 
x-coordinate  is  normal  to  the  absorbing  surface.  In  2-D,  the  paraxial  equation  is, 

Ett±cEtx-^(Eyy)  =  0,  (28) 

where  c  =  1  /y'ep  denotes  the  wavespeed,  while  in  3-D  this  becomes 

Ett±cEtx-£(Eyy+Ea)  =  0.  (29) 

The  purpose  of  the  exercise  is  to  derive  a  consistent  boundary  condition  from  these  paraxial 
approximations  of  Maxwell’s  equations.  This  is  accomplished  by  applying  the  standard  Galerkin  finite 
element  formulation  to  equations  (28-29),  from  which  an  implementation  of  the  boundary  conditions 
consistent  with  the  interior  domain  may  be  derived.  However,  there  does  not  appear  to  be  a  single, 
consistent  way  of  doing  this,  particularly  in  the  context  of  the  reduced  integration  techniques  used 
in  finite  element  algorithms. 

For  the  2-D  case,  no  difficulties  are  encountered.  Clearly,  more  work  on  the  implementation  of 
paraxial  boundary  conditions  in  finite  element  algorithms  is  required.  It  appears  that  a  complete 
paraxial  finite  element  formulation  with  a  layer  of  paraxial  elements  around  the  model  instead  of 
merely  paraxial  boundaries  will  be  a  useful  direction  for  development. 
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Although  effective,  the  4th  order  paraxial  conditions  are  not  naturally  compatible  with  discrete 
algorithms,  nor  do  they  readily  admit  conductivity.  The  problem  is  that  these  methods  are  based  on 
an  analytical  approximation  of  Maxwell’s  equations  rather  than  on  an  approximation  of  the  discrete 
finite  element  form  of  these  equations.  There  is  a  subtle  but  real  difference.  Approaches  that  attack 
the  discrete  equations  themselves  have  been  studied  recently.  In  particular,  Sandler  (1991)  developed 
a  mechanically  based  absorbing  condition  for  nonlinear  solid  mechanics  applications.  This  condition 
is  superior  to  the  paraxial  approach  described  previously,  because  it  eliminates  wavespeed  from  the 
formulation  (which  in  nonlinear  calculations  is  problematic)  and  provides  better  directionality  than 
the  2nd  order  paraxial  or  normal  incidence  condition. 

9.  Time  Domain  vs.  Frequency  Domain 

Time-domain  methods  currently  provide  the  quickest,  least  computer  memory  intensive,  most 
robust  path  to  a  solution.  Advances  in  numerical  methods  may  someday  change  this  answer.  There 
do  exist  reasons  for  using  frequency-domain  simulations  if  they  become  competitive  in  CPU  and 
memory  cost. 

The  primary  advantage  of  the  time-domain  solver  is  that  it  embodies  an  efficient  explicit  algorithm 
that  requires  minimal  memory,  thus  permitting  solutions  of  the  largest  finite  element  model  possible 
on  any  given  machine.  It  is  also  robust  and  deterministic  in  the  sense  that  if  run  long  enough, 
steady-state  will  be  achieved  and  a  solution  will  be  found.  Additional  advantages  of  the  time-domain 
approach  are  that  it  is  directly  extensible  to  nonlinear  problems  where  material  properties  change  with 
time,  e.g.,  the  bleaching  process  in  photolithography,  wave  arrivals  may  be  separated  in  time 
providing  better  insight  into  physical  processes,  and  transient  (pulsed)  or  other  nonsinusoidal  signals 
are  easily  accommodated. 

Disadvantages  of  time-domain  solvers  are  that  user  intervention  is  required  to  assess  when 
steady-state  has  been  reached.  One  can  envision  problems  where  a  long  simulation  may  be  necessary 
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to  achieve  steady-state.  Also,  the  steady-state  quantities  (amplitude  and  phase)  are  not  primary  in 
the  time  domain  and  must  be  obtained  by  a  secondary  calculation  after  steady- state  has  been  achieved. 

Given  that  we  are  looking  for  a  frequency-domain  solution,  the  advantages  of  a  straightforward 
frequency-domain  formulation  appear  obvious.  The  difficulty  is  that  for  realistically  sized  models  this 
formulation  requires  the  solution  of  an  enormous,  sparse,  linear  system  of  equations.  In  addition,  the 
linear  system  has  undesirable  numerical  properties,  namely,  it  is  complex,  non-Hermitian  (due  to 
absorbing  boundary  conditions  and/or  conductivity),  nonsymmetric  (due  to  absorbing  boundary 
conditions),  and  is  typically  indefinite. 

There  are  two  basic  methods  of  solving  linear  systems:  direct  (some  variation  of  Gaussian 
elimination)  and  iterative.  Direct  methods  are  typically  divided  into  two  phases:  factorization  of  the 
coefficient  matrix,  followed  by  back-substitution  to  obtain  the  solution.  The  payoff  here  is  that  most 
of  the  work  is  performed  in  the  factorization,  so  that  solutions  for  additional  right-hand  sides,  e.g., 
different  illumination  sources,  can  be  obtained  from  a  relatively  cheap  backsolve.  The  drawback  of 
direct  methods  is  that  they  require  large  amounts  of  memory  and  CPU  time.  For  example,  consider 
a  100  x  100  x  100  element  3-D  model  with  three  electric  field  components  unknown  at  each  node, 
in  complex  arithmetic.  This  is  a  modest  model  size.  For  this  problem,  a  standard  band  solver  such 
as  ZBGFA  out  of  UNPACK  (Dongarra  et  aL  1979)  requires  5.73  x  10n  words  of  storage,  which  far 
exceed  the  capacity  of  any  present  day  computer.  The  number  of  floating  point  operations  to  solve 
such  a  system  is  also  prohibitive.  It  is  worth  noting  here  that  this  linear  system  is  very  sparse,  i.e., 
only  2.5  x  108  or  0.04%  of  the  entries  are  nonzero.  Observe  that  even  the  nonzero  entries  exceed 
the  largest  available  machines,  but  only  by  a  small  margin.  Some  work  has  been  done  on  sparse 
system  direct  solvers  that  attempt  to  economize  on  storage  relative  to  the  basic  band  solvers.  Thus, 
we  could  conservatively  concur  with  the  conventional  wisdom  which  holds  that  3-D  problems  cannot 
be  solved  directly. 

The  alternative  to  direct  solvers  is  iterative  methods  in  which  an  initial  guess  for  the  solution  (zero 
if  nothing  better  is  available)  is  successively  refined  until  the  error  becomes  “small.”  Conventional 
wisdom  says  that  the  Krylov  subspace-type  methods  are  best,  which  includes  the  conjugate  gradient 
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approach.  Frequently,  some  type  of  preconditioner  is  used  in  conjunction  with  these  methods  to 
accelerate  convergence.  For  some  classes  of  matrices,  e.g.,  positive  definite  symmetric,  these 
methods  work  extremely  welL  The  major  portion  of  the  CPU  effort  is  in  computing  a  product  of  the 
coefficient  matrix  with  a  vector  and  in  computing  inner  products  of  two  vectors.  Thus,  the 
insurmountable  memory  requirements  of  the  direct  methods  can  be  avoided.  Additional  vectors  of 
length  N  (the  total  number  of  unknowns),  which  form  the  basis  of  the  Krylov  Subspace,  are  usually 
required. 

Things  to  consider  in  choosing  an  iterative  method  are:  1)  the  total  amount  of  work  required  by 
the  method  and  the  preconditioner  to  achieve  an  acceptable  solution;  2)  the  total  amount  of  memory 
required  by  the  method  and  the  preconditioner;  3)  the  possibility  of  breakdown-the  algorithm  foils, 
or  hangs-the  convergence  rate  becomes  so  low  that  a  solution  is  not  attainable;  and  4)  determinism-a 
solution  might  be  obtainable  only  by  trial-and-error  tweaking  of  the  preconditioner  or  iterative  solver 
parameters. 

10.  Summary 

Recently,  Freund  (1991 , 1992)  proposed  a  new  Krylov-subspace  method  that  he  calls  QMR  (for 
quasi-minimal  residual)  and  extended  it  to  general  nonsingular,  non-Hermitian  systems.  He 
specifically  targets  the  systems  arising  from  frequency-domain  discrete  numerical  methods.  He  has 
also  devised  a  strategy  that  circumvents  one  cause  of  breakdown,  though  incurable  breakdowns  are 
still  possible  in  theory. 

For  frequency-domain  finite  element  solvers,  solution  techniques  for  the  resulting  class  of  linear 
systems  are  still  in  the  research  phase.  Significant  advances  have  been  achieved  in  recent  years,  and 
more  advances  are  expected.  Trying  new  solution  algorithms  is  a  worthwhile  exercise,  which 
provides  feedback  to  the  mathematicians  and  computer  scientists,  and  is  to  everyone’s  advantage. 
On  the  other  hand,  it  is  premature  to  promote  such  solvers  as  production  level  tools  for  engineers. 
If  large-scale  calculations  need  to  be  done  today,  time-domain  techniques  provide  the  most  practical 
means  of  doing  them. 
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